clear all
close all

% =====================================================================
% ====  "Risk-Sharing in Village Economies Revisited", by Bold and Broer

% ====   This programme solves the limited commitment village economy for
% both the standard specification (individual deviations) and the
% alternative specification (coalitional deviations)

% This version considers preference heterogeneity

% ====    Written by Tessa Bold and Tobias Broer

% ====    This version Feb 2021
% =====================================================================


filetosave='Heterog_Vilcoal_20_Oct_vill_cd_incproc_outl'; % name of workspace saved

%====================================================
% select parameters
% ===================================================

% ======== a. individual's income process ========

village=3;                         % village indicator: 1 - Aurepalle; 2-Kanzara; 3-Shirapur
coaldev=0;                         % set to 1 for coalitional deviations
incproc=1;                        % 1 to 4 income processes define size of measurement error in income
heterog=1;                         % set to 1 if you want heterogeneous risk-aversion
buildup=0;    %set to 1 and change Qrho and Qbeta loop to calculate for all group sizes
if heterog==1
    coaldev=0;  buildup=0;
end
unconditional=0;                         % only applies to village 1: set this to 0 to load datasheet without outlier correction, 1 with larges observations eliminated, 2 with first period eliminated altogether

filetosaveest=['Momentsest' int2str(village) '410'];
Tessa1Tobi2=2;

server=0;				% set to 1 to run on server; o wise path for workstations is used
if village==1
    vss_high_inddev=34-1;                 % the number of total village members is vss_high_inddev+1 under individual deviations
    vss_high_coaldev=34-1;
elseif village==2;
    vss_high_inddev=37-1;                 % the number of total village members is vss_high_inddev+1 under individual deviations
    vss_high_coaldev=37-1;
elseif village==3
    vss_high_inddev=31-1;                 % the number of total village members is vss_high_inddev+1 under individual deviations
    vss_high_coaldev=31-1;
end
vssbuildupvec_inddev=[1:9 12 17 24 vss_high_inddev]; % if you choose buildup==1, rhis is the vector of village sizes for which the solution is calculated in the individual deviations case
if server==1
    outputpath=sprintf('/home/workgroups/testob/Output estimation June 2018 conditional');
    programpath=sprintf('/home/workgroups/testob');
else
    if Tessa1Tobi2==2
        datapath='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Data\Output';
        outputpath=('C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Programmes\Main Results and Robustness')
        programpath=('C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Programmes')
    elseif Tessa1Tobi2==1
        
        outputpath=sprintf('C:/Users/tbold/Dropbox/Tessa and Tobi/Programmes')
        programpath=sprintf('C:/Users/tbold/Dropbox/Tessa and Tobi/Programmes')
        
    end
end
betavec=[[0.5:0.04:0.7],[0.72:0.02:0.98]];
rhovec=[0.5,1.5];
rho_rov=1;

%main execution options
gapcritset=0.0005%0.0002; %convergence criterion for change in values

% other execution options
Nincproc=1;                         % number of income processes to look at
T=6200; Tinit=201;                 % number of simulation periods in total, and number of initial periods to discard
maxcerr=0.9;                          % max cons meas error in log levels
N_cerr=30;                           % number of support points for cons meas error
maxincerr=2/3;                      % max inc meas error in multiples of inc variance (in levels)
% set to one if you want to specify the rest of village in terms of average
% consumption
rov_average=1;
% set this to 1 if you want to calculate the moments for levels of
% consumption based on log consumption
momentclevinlog=1;
Rov_sb=2; % set this to 1 if you want the outside option for the Rov in the case of coalitional deviations to equal the average value from the solution to the next smaller village (second-best), rather than
% just the utility from consuming average ROV income (first-best)
% convergence criterion and maximum number of iterations
maxitgapset=400;
% this is a punishment factor: default involves a consumption penalty in the first period equal
% to caut*Pun_fac
Pun_fac=0.0;
% grid of time-varying planner weights
%lambda=[0.1,0.125,0.15,0.175,0.2,0.25,0.3,0.35:0.05:0.65,0.7,0.75,0.8,0.825,0.85,0.875,0.9]';
lambda=[0.05:0.001:0.95]';%linspace(0.05,0.95,301)';
m=length(lambda);
if floor((m+1)/2)~=((m+1)/2)
    disp('care: lambda has to have uneven no of points in first dimension')
    stop
end
S_KEEP=nan(1000*3,2,vss_high_inddev);                  % pre-allocate the matrices for coalitional deviations, so that you don't run simulations vss times
S_rov_KEEP=nan(100,1,vss_high_inddev);              % this one is for a maximum of vss=19 for the CD
P_rov_KEEP=nan(100,1000,vss_high_inddev);
PHI_KEEP=nan(m,1000*3,2,vss_high_inddev);
IDIO_DIST_ROV_KEEP=nan(1000*3,3,vss_high_inddev);

% ====================================================
% I. Specify Income process for individual and Rest of Village (rov)
% ====================================================
% data moments for log income from ICRISAT data - all villages
% first col: Aurepalle, secdon: Kanzara, third: Shirapur
if unconditional==0
    datafile='armoments_fe0';
elseif unconditional==1
    datafile='armoments_fe1';
end
eval(['load ' datafile '.out;']);
eval(['covyylag=' datafile '(1,:);']);
eval(['vardy=' datafile '(2,:);']);
eval(['vary=' datafile '(3,:);']);

if Nincproc>1
    varnu_max=maxincerr*vardy/2;                % nu is measurement error in log levels
else
    varnu_max=zeros(1,3);
end
rho_min=covyylag(village)/vary(village);                      % AR(1) coefficient w/o meas error
rho_max=covyylag(village)/(vary(village)-varnu_max(village));          % AR(1) coefficient w max meas error
rho1vec=linspace(rho_min,rho_max,Nincproc);
varnu=vary(village)-covyylag(village)./rho1vec;
vareps=(vary(village)-varnu).*(1-rho1vec.^2);

% ===================================
% iterate over income processes
% ===================================
disp('now solving for village,coaldev,incproc, unconditional - press any key to continue')
disp([village coaldev unconditional])
pause

%Momentsest=nan(28,11,size(betavec,2),size(rhovec,2),vss_high_inddev);
indicator1=zeros(1,size(rhovec,2));
%for incproc=1:3
% specify a markov process for individual income
[incomevals, Q1] = rouwenhorst(rho1vec(incproc),vareps(incproc)^0.5,3);
incomevals=exp(incomevals');
cumQ=cumsum(Q1,2);

% set minimum and maximum village size
if coaldev==1                       % for caolitional deviations,
    vss_set  =[1:19 19+round((vss_high_coaldev-19)/2) vss_high_coaldev];
    vss_high=vss_set(end);
else
    if buildup==0% for individual deviations, start with small village and build up, but coarser
        if heterog~=1
            vss_set=[vss_high_inddev];
            lagind=1;
            vss_high=vss_set(end);% an indicator
        else
            % for heterogeneity, number of individuals must be multiple of
            % number of groups
            vss_high_inddev=ceil(vss_high_inddev/length(rhovec))*length(rhovec);
            vss_set=vss_high_inddev-1;
            lagind=1;
            vss_high=vss_set(end);% an indicator
        end
    elseif buildup==1
        vss_set=vssbuildupvec_inddev;
        lagind=1;
        vss_high=vss_set(end);% an indicator
    end
end

Income0=nan(vss_high_inddev+1,T,size(betavec,2),size(rhovec,2),Nincproc);
Csim0=nan(vss_high_inddev+1,T,size(betavec,2),size(rhovec,2),Nincproc);

% ===================================
% iterate over discount factors
% ===================================
for Qbeta=1:size(betavec,2) %[1,3,5,7,9,11,13,14,2,4,6,8,10,12]%
    
    % ===================================
    % iterate over CRRA
    % ===================================
    waut_lag=[]; waut_avg_lag=[];
    for vsscount=1:length(vss_set)       % iterate over village size
        
        vss=vss_set(vsscount);
        vss_max=max(vss_set);
        if vsscount>1 & coaldev==1
            lagind=lagind+(vss-(vss_set(vsscount-1)+1));
        end
        
        disp('now solving for parameters')
        disp([Qbeta])
        disp('vss')
        disp(vss)
        % this is the scaling factor: if you want RoV consumption as average consumption per HH, scale rov variables
        % by vss
        if rov_average==1
            ScaleRov=vss;
        else
            ScaleRov=1;
        end
        csim=[];  cc =[]; eeewaut =[]; eewaut =[]; ewaut =[]; waut =[]; rr =[]; cshare=[]; raw =[]; cshare =[]; dcshare =[]; dcshare_raw =[]; S_rov =[]; P_rov =[];
        S_rov_new =[]; P_rovnew =[]; S =[]; P =[]; caut =[]; c =[]; waut =[]; gammagrid =[]; gammar =[]; phi =[]; phisol =[]; csol =[];csolold=[]; cfirstbest =[];
        wfirstbest =[]; ewfirstbest =[]; wsbnewsol=[];
        Y =[]; vs=[];  idio_dist_autnew=[];idio_dist_aut_new=[];normperf=[];w1=[];indSj=[]; indSjj=[]; Ssim=[]; aggvillincome=[];aggincgrowth=[];aggvillincgrowth=[];
        S_rov_aut=[]; S_rov_sb=[];
        w =[];  phinew =[]; Yagg =[]; wint =[]; wsb =[]; ewsb =[]; csol =[]; phisol =[]; nonconv =[];   gammarfun =[];  x0 =[]; index =[]; ewsb =[]; bind=[];
        
        
        % ===================================
        % build the state space of ROV income
        % ===================================
        idio_dist=[];
        idio_dist_aut=[];
        idio_dist_sb=[];
        
        for k1=0:vss
            for k2=0:vss
                for k3=0:vss
                    II=zeros(1,3);
                    if k1+k2+k3==vss
                        % the vector of aggregate incomes from
                        % permutations of vss households
                        S_rov=[S_rov; k1*incomevals(1,1) + k2*incomevals(2,1) + k3*incomevals(3,1)];
                        % the corresponding distribution of
                        % individual income values
                        idio_dist=[idio_dist; k1 k2 k3];
                        
                        % this computes the 'best' distribution of incomes
                        % of size vss-lagind
                        if vss>1
                            II=zeros(1,3);
                            temp=[k1 k2 k3];
                            % this takes away the lowest incomes to
                            % derive the distribution of incomes of the
                            % 'best' coalition lagind individuals
                            % This coalition is one smaller than the
                            % largest sustainable coalition because it
                            % will consist of 1 agent and m-1 guys from
                            % the largest sustainable coalition of size
                            % m
                            
                            for jj=1:lagind
                                II(min(find(temp>0)))=II(min(find(temp>0)))+1;
                                temp=[k1 k2 k3]-II;
                            end
                            idio_dist_aut=[idio_dist_aut; [k1 k2 k3]-II];
                            S_rov_aut=[S_rov_aut; ([k1 k2 k3]-II)*incomevals];
                            II=zeros(1,3);
                            temp=[k1 k2 k3];
                            % this takes away the lowest incomes to
                            % derive the distribution of incomes of the
                            % 'best' coalition
                            % consisting of m individuals from only the
                            % rest of village
                            for jj=1:lagind-1
                                II(min(find(temp>0)))=II(min(find(temp>0)))+1;
                                temp=[k1 k2 k3]-II;
                            end
                            idio_dist_sb=[idio_dist_sb; [k1 k2 k3]-II];
                            S_rov_sb=[S_rov_sb; ([k1 k2 k3]-II)*incomevals];
                        end
                    end
                end
            end
        end
        
        % resort S_rov in ascending order
        [S_rov,S_rovind]=sort(S_rov);
        idio_dist=idio_dist(S_rovind,:);
        
        if vss>1
            idio_dist_aut=idio_dist_aut(S_rovind,:);
            S_rov_aut=S_rov_aut(S_rovind);
            idio_dist_sb=idio_dist_sb(S_rovind,:);
            S_rov_sb=S_rov_sb(S_rovind);
        end
        
        
        % ===================================
        % simulate aggregate income transitions
        % ===================================
        tic
        T_trans=50000;
        aggincind=zeros(length(S_rov),length(S_rov));
        for jjj=1:length(S_rov)
            jjj;
            inccount=zeros(T_trans,vss,3);
            agginc=zeros(T_trans,3);
            for iii=find(idio_dist(jjj,:)>0)
                
                for kkk=1:idio_dist(jjj,iii)
                    shock=rand(T_trans,1);
                    rr=(find(shock<=cumQ(iii,1)));
                    inccount(rr,kkk,iii)=1;
                    rr=(find(cumQ(iii,1)<shock & cumQ(iii,2)>=shock));
                    inccount(rr,kkk,iii)=2;
                    rr=(find(cumQ(iii,2)<shock & cumQ(iii,3)>=shock));
                    inccount(rr,kkk,iii)=3;
                end
                
                agginc(:,iii)=sum(incomevals(inccount(:,1:idio_dist(jjj,iii),iii)),2);
            end
            agginc=sum(agginc,2);
            for jjj2=1:length(S_rov)
                qqq=(agginc<=S_rov(jjj2)+0.0000000001).*(agginc>=S_rov(jjj2)-0.0000000001);
                aggincind(jjj,jjj2)=sum(qqq);
            end
            P_rov(jjj,:)=aggincind(jjj,:)/sum(aggincind(jjj,:));
        end
        disp('Simulation for agg income process takes')
        toc
        
        
        %scale everything correctly
        if vss>1
            S_rov_aut=S_rov_aut/(max(ScaleRov-lagind,1));
            S_rov_sb=S_rov_sb/(max(ScaleRov-lagind+1,1));
        end
        
        S_rov=S_rov/ScaleRov;
        
        % ===================================
        % state space: combinations of indiv and village income
        % ===================================
        S=[kron(incomevals,ones(size(S_rov))),kron(ones(size(incomevals)),S_rov)];
        S_rov_aut=[kron(ones(size(incomevals)),S_rov_aut)];
        S_rov_sb=[kron(ones(size(incomevals)),S_rov_sb)];
        P=kron(Q1,P_rov);
        
        Y=S(:,1)+ScaleRov*S(:,2);				% possible aggregate income outcome
        state=length(S);                    % number of states
        agent=2;                            % number of 'agents' in the HH vs. RoV approximation equals 2
        STATE(vss)=state;
        for Qrho=1:size(rhovec,2)
            
            beta=betavec(Qbeta)
            rho=rhovec(Qrho)
            
            if coaldev==1
                lagind=[];
            else
                lagind=1;
            end
            % ===================================
            % Compute set of autarky values
            % ===================================
            % autarky values: coalitional deviations imply that the individual's
            % outside option is given by utility of caut for 1 period, and then the
            % expected value being in the 'best' insurance coalition of
            % size vss-1, with most 'equal' weights, i.e.
            % gamma=1/(vss-1) or, if this does not yield the autarky value corresponding to a coalition of vss-2, the value that is closest to it
            for k=1:agent
                caut(:,k)=S(:,k);
            end
            waut=zeros(state,agent);
            waut_avg_lag=zeros(state,agent);
            waut_sb=zeros(state,agent);
            % this is just the PDV of utility from consumption of
            % income caut
            eeewaut(:,1)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,1));
            eeewaut(:,2)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho_rov,caut(:,2));
            
            % in first iteration, or for individual deviations, waut is just consumption of income
            if vss==1 || coaldev==0
                for j=1:state
                    waut(j,1)=utility(rho,caut(j,1)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,1);
                    % autarky values for the rest of village are always the same - perfect
                    % risk-sharing among the rest of village
                    waut(j,2)=utility(rho_rov,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,2);
                end
                
            else
                % this calculates autarky values for the coalitional
                % deviations case
                %if you were going to use the second-best constrained
                %value, you first need to create the average lagged
                %utility weighted according to income of agent 1 and rest
                %of the village from the previous iteration
                %1) look for the entry in S_rov_sb, which has the same
                %average/aggregate income (1*S_lag(:,1)+(vss-lagind)*S_lag(:,2))
                %2) this is usually more than one entry, since S_lag
                %differs by which income agent 1 has and which income
                %is rest of village. Here, we don't care about that
                %anymore, only about the aggregate. Therefore we also
                %need to average utility (hardly matters for which
                %point you do this, to be really precise should take
                %half form both), but important is to average utility
                %as (waut_lag(ind_sb,1)+(vss-lagind)*waut_lag(ind_sb,2))/(vss-lagind+1)
                %NOTE THIS IS WRITTEN FOR AVERAGE INCOME IN REST OF
                %VILLAGE; NOT AGGREGATE, THE AVERAGING WOULD NEED
                %%CHANGING THEN.
                
                for j=1:state
                    % this indicates combinations of individual and
                    % rest-of-group incomes in the next smaller stable
                    % coalition that have the same total income as S_rov_sb(j,2)
                    S_lag_ind_sb=find((S_rov_sb(j,1)<=(S_lag(:,1)+S_lag(:,2)*(vss-lagind))/(vss-lagind+1)+0.00001 & S_rov_sb(j,1)>=(S_lag(:,1)+S_lag(:,2)*(vss-lagind))/(vss-lagind+1)-0.00001));
                    % this calculates the autarky utility for the
                    % individual and the rest-of-group for those
                    % combinations - one period of autarky plus expected
                    % disc utility from tomorrow
                    wwww=[utility(rho,S_lag(S_lag_ind_sb,1))+beta*P_lag(S_lag_ind_sb,:)*waut_lag(:,1) utility(rho_rov,S_lag(S_lag_ind_sb,2))+beta*P_lag(S_lag_ind_sb,:)*waut_lag(:,2)];
                    % deviation utility of rest-of-group today is the
                    % average of individual and rest-of-group utilities for
                    % next smaller stable group
                    % equal-prob draw across those states implies you just
                    % take the mean of the mean
                    wwwalt(j)=mean(mean(wwww,2));
                    www=(waut_lag(S_lag_ind_sb,1)+(vss-lagind)*waut_lag(S_lag_ind_sb,2))/(vss-lagind+1);
                    %%taking the average here is not quite correct,
                    %%clearly this should be weighted by how likely either
                    %%state is to occur given the previous state, but that
                    %%just seems needlessly complicated.
                    waut_avg_lag(j,2)=mean(www);
                    %waut_avg_lag(j,2)=wwwalt(j); %Tobi's version, makes little difference, all seems fine
                    
                end
                
                
                for j=1:state
                    % indicator for the 'autarky' state in the state
                    % space of the previous iteration S_lag
                    S_lag_2=(S_rov_aut(j)==S_lag(:,2)) ;
                    S_lag_1=(S(j,1)==S_lag(:,1));
                    S_lag_ind=find(S_lag_1.*S_lag_2==1);
                    % for the individual, autarky value is utility of
                    % current consumption plus relevant trans
                    % probability times waut_lag, the vector of values
                    % from the next smaller insurance group
                    waut(j,1)=utility(rho,caut(j,1)*(1-Pun_fac))+beta*P_lag(S_lag_ind,:)*waut_lag(:,1);
                    % for rest of village, the autarky value is just
                    % expected disc utility of current consumption
                    
                    waut(j,2)=utility(rho_rov,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,2);
                    
                    %%and the second best utilities
                    waut_sb(j,1)=waut(j,1);
                    waut_sb(j,2)=utility(rho_rov,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*waut_avg_lag(:,2);
                    %waut_sb_2(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P_rov(j,:)*waut_avg_lag(1:6,2);
                    %JUST A SANITY CHECK
                end
                
            end
            
            if Rov_sb==1 & vss>1 & coaldev==1
                waut=waut_sb;
            end
            
            
            % ====================================================
            % Compute constrained insurance
            % ====================================================
            gammagrid=[lambda vss/ScaleRov*(1-lambda)];
            gammagrid=[gammagrid ones(m,1) gammagrid(:,2)./gammagrid(:,1)];
            gammagridold=gammagrid;
            clear gammagrid
            gammagrid=[linspace(1.2*max(S(:,1)),min(S(:,1)),m)',linspace(min(S(:,2)),1.2*max(S(:,2)),m)',ones(m,1),linspace(1.2*max(S(:,1)),min(S(:,1)),m)'./linspace(min(S(:,2)),1.2*max(S(:,2)),m)'];
            
            % ====== make grids of weights, consumption, and Lagrange multipliers for updating ======
            n=length(gammagrid);
            w=zeros(n,state,agent);
            gammar=zeros(n,state,agent);
            c=zeros(n,state,agent);
            phi=zeros(n,state,agent);
            
            % 'first-best' consumption given aggregate resources and planner weights
            for j=1:state
                c(:,j,1)=Y(j)./(1+ScaleRov*gammagrid(:,4));
                c(:,j,2)=(Y(j)-c(:,j,1))/ScaleRov;
            end
            
            %for j=1:state
            % c(:,j,1)=Y(j)*2/(vss+1)*gammagrid(:,1);
            % c(:,j,2)=Y(j)*2/(vss+1)*gammagrid(:,2);
            %end
            
            cfirstbest=c;
            %  relative planner weights
            for j=1:state
                for k=1:agent
                    gammar(:,j,k)=gammagrid(:,agent+k);
                end
            end
            
            % Compute the first best value functions without binding
            % constraints - i.e. cfirstbest for ever
            w=zeros(n,state,agent);
            
            dww=1;
            wold=1/(1-beta)*utility(rho,cfirstbest(:,:,:));
            while dww>0.00001
                w(:,:,1)=utility(rho,cfirstbest(:,:,1))+beta*wold(:,:,1)*P';
                w(:,:,2)=utility(rho_rov,cfirstbest(:,:,2))+beta*wold(:,:,2)*P';
                
                dww=max(max(max(abs(w-wold))));
                wold=w;
            end
            wfirstbest=w;
            
            % continuation utility as function of planner weight and state of the
            % economy
            ewfirstbest(:,:,1)=wfirstbest(:,:,1)*P';
            ewfirstbest(:,:,2)=wfirstbest(:,:,2)*P';
            
            %disp('384')
            % ====================================================
            % Start the loop with the enforcement constraints
            % ====================================================n
            ewsb=ewfirstbest;
            wsb=wfirstbest;
            options=optimset('Display','off');
            eps=0.0000001;t=0;itgap=0; gap=2; gaphist=gap; gapincrease=0;
            
            csolold=cfirstbest;
            if beta<0.85
                gapcrit=gapcritset;
                maxitgap=maxitgapset;
            else
                gapcrit=gapcritset*1
                maxitgap=maxitgapset;
            end
            
            while gap >=gapcrit && itgap<=maxitgap
                rr=[];
                rr1=[];
                rr2=[];
                rrnone=[];
                rrboth=[];
                nonconv=zeros(state,3)*nan;
                genuine=zeros(state,1)*nan;
                itgap=itgap+1;
                coalnonsust=nan(state,2);
                count=0; countgrid=[ones(n,state)];
                for j=1:state
                    %instead of solving problem for each grid point,
                    %note that there are only three possibilities: 1)
                    %Constraint binding for Person 1 only, the
                    %solution is the same for each lambda for which
                    %this is the case
                    index=[1 j];
                    
                    rr1=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)<waut(j,1) & utility(rho_rov,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2)) ;
                    
                    if size(rr1)>0
                        bind=1;
                        rr=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)>=waut(j,1));
                        if isempty(rr)==1
                            coalnonsust(j,1)=-2;
                            output=[nan,nan,-50];
                        else
                            clear x0
                            
                            x0(1) = cfirstbest(rr(1),j,1);  %redundant because we no longer need an initial guess
                            %T vectors of consumption and utilty that meet
                            %participation constraints
                            [x, output]=agent1maxrho_rov_het_06_Mar(x0,index,ewsb,waut,bind,Y,n,beta,rho,rho_rov,ScaleRov,cfirstbest,rr);
                            
                            %T 29_Mar This keeps track of nonconvergence, but
                            %allows the programme to continue - sometimes it
                            %converges in a subsequent iteration
                        end
                        if output(3) <=0
                            genuine(j,1)=(rr1(end)-rr(1));
                            nonconv(j,1)=output(3);
                            csol(rr1,j,1)=ones(size(rr1,1),1)*caut(j,1);
                            csol(rr1,j,2)=ones(size(rr1,1),1)*caut(j,2);
                            gammar(rr1,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                            wsbnewsol(rr1,j,1)=ones(size(rr1,1),1)*waut(j,1);
                            wsbnewsol(rr1,j,2)=ones(size(rr1,1),1)*waut(j,2);
                            %phisol(rr1,j,1)=gammagrid(rr1,4).^rho./gammar(rr1,j,2).^rho-1;
                            phisol(rr1,j,1)=(cfirstbest(rr1,j,2).^rho_rov./cfirstbest(rr1,j,1).^rho)./(csol(rr1,j,2).^rho_rov./csol(rr1,j,1).^rho)-1;
                            phisol(rr1,j,2:agent)=0	;
                            
                            
                        else
                            
                            
                            csol(rr1,j,1)=x(1);
                            gammar(rr1,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                            csol(rr1,j,2)=(Y(j)-x(1))/ScaleRov;
                            wsbnewsol(rr1,j,1)=output(1);
                            wsbnewsol(rr1,j,2)=output(2);
                            %phisol(rr1,j,1)=gammagrid(rr1,4).^rho./gammar(rr1,j,2).^rho-1;
                            phisol(rr1,j,1)=(cfirstbest(rr1,j,2).^rho_rov./cfirstbest(rr1,j,1).^rho)./(csol(rr1,j,2).^rho_rov./csol(rr1,j,1).^rho)-1;
                            phisol(rr1,j,2:agent)=0	;
                        end
                    end
                    
                    
                    
                    % 2. Constraint binding for agent 2 (RoV) only
                    rr2=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)>=waut(j,1) & utility(rho_rov,cfirstbest(:,j,2))+beta*ewsb(:,j,2)<waut(j,2)) ;
                    
                    if size(rr2)>0
                        bind=2; %i.e. agent 2 has the binding constraint
                        rr=find(utility(rho_rov,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2));
                        if isempty(rr)==1
                            coalnonsust(j,2)=-2;
                            output=[nan,nan,-50];
                        else
                            clear x0
                            x0(1) = cfirstbest(rr(end),j,1); %redundant because we no longer need an initial guess
                            
                            [x, output]=agent1maxrho_rov_het_06_Mar(x0,index,ewsb,waut,bind,Y,n,beta,rho,rho_rov,ScaleRov,cfirstbest,rr);
                        end
                        % output(3) equals 1 if constraints
                        % are satisfied, but -2 if not
                        if output(3)<=0
                            nonconv(j,2)=output(3)*100;
                            csol(rr2,j,1)=ones(size(size(rr2,1),1),1)*caut(j,1);
                            gammar(rr2,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                            csol(rr2,j,2)=ones(size(size(rr2,1),1),1)*caut(j,2);
                            wsbnewsol(rr2,j,1)=ones(size(rr2))*waut(j,1);
                            wsbnewsol(rr2,j,2)=ones(size(rr2))*waut(j,2);
                            %phisol(rr2,j,2)=gammar(rr2,j,2).^rho_rov./gammagrid(rr2,4).^rho_rov-1;
                            phisol(rr2,j,2)=(csol(rr2,j,2).^rho_rov./csol(rr2,j,1).^rho)./(cfirstbest(rr2,j,2).^rho_rov./cfirstbest(rr2,j,1).^rho)-1;
                            phisol(rr2,j,1)=0;
                        else
                            
                            
                            csol(rr2,j,1)=x(1);
                            gammar(rr2,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                            csol(rr2,j,2)=(Y(j)-x(1))/ScaleRov;
                            wsbnewsol(rr2,j,1)=output(1);
                            wsbnewsol(rr2,j,2)=output(2);
                            %phisol(rr2,j,2)=gammar(rr2,j,2).^rho_rov./gammagrid(rr2,4).^rho_rov-1;
                            phisol(rr2,j,2)=(csol(rr2,j,2).^rho_rov./csol(rr2,j,1).^rho)./(cfirstbest(rr2,j,2).^rho_rov./cfirstbest(rr2,j,1).^rho)-1;
                            phisol(rr2,j,1)=0;
                        end
                        
                    end
                    
                    % 3. Constraint all slack
                    rrnone=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)>=waut(j,1) & utility(rho_rov,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2)) ;
                    
                    if size(rrnone)>0
                        
                        csol(rrnone,j,1:agent)=cfirstbest(rrnone,j,1:agent);
                        gammar(rrnone,j,1:agent)=gammagrid(rrnone,agent+1:2*agent);
                        phisol(rrnone,j,1:agent)=0;
                        
                        %for g=1:agent
                        %   wint(g)=interp1(gammagrid(:,4),ewsb(:,j,g),gammagrid(i,4));
                        %end;
                        wsbnewsol(rrnone,j,1)=utility(rho,cfirstbest(rrnone,j,1))+beta*ewsb(rrnone,j,1);
                        wsbnewsol(rrnone,j,2)=utility(rho_rov,cfirstbest(rrnone,j,2))+beta*ewsb(rrnone,j,2);
                        
                        
                        
                    end
                    
                    % 3. Constraint all binding, this can happen for
                    % coalitional deviations
                    rrboth=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)<waut(j,1) & utility(rho_rov,cfirstbest(:,j,2))+beta*ewsb(:,j,2)<waut(j,2)) ;
                    
                    if size(rrboth)>0 & itgap==1
                        
                        nonconv(j,3)=-20000;
                        csol(rrboth,j,1)=caut(j,1);
                        gammar(rrboth,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                        csol(rrboth,j,2)=caut(j,2);
                        wsbnewsol(:,j,1)=waut(j,1);
                        wsbnewsol(:,j,2)=waut(j,2);
                        %%note I am not fixing this now for the
                        %%heterogeneity, should not happen for ID.
                        %phisol(rrboth,j,2)=gammar(rrboth,j,2).^rho./gammagrid(rrboth,4).^rho_rov-1;
                        %phisol(rrboth,j,1)=0;
                    end
                    
                    
                end
                
                gapold=max(norm(wsbnewsol(:,:,1)-wsb(:,:,1)),norm(wsbnewsol(:,:,2)-wsb(:,:,2)));
                %gap=max(max(max(abs(wsbnewsol(:,:,1)-wsb(:,:,1)))),max(max(abs(wsbnewsol(:,:,2)-wsb(:,:,2)))))
                gap1=max(max(max(abs(invutility(beta,rho,wsbnewsol(:,:,1))-invutility(beta,rho,wsb(:,:,1)))./invutility(beta,rho,wsb(:,:,1)))),max(max(abs(invutility(beta,rho_rov,wsbnewsol(:,:,2))-invutility(beta,rho_rov,wsb(:,:,2)))./invutility(beta,rho_rov,wsb(:,:,2)))))/(1-beta);
                gap2=max(max(max(abs(csol(:,:,1)-csolold(:,:,1))./csol(:,:,1))),max(max(abs(csol(:,:,2)-csolold(:,:,2))./csol(:,:,2))));
                gap=max(gap1,gap2);
                gaphist(itgap+1)=gap;
                disp('gap,gapold')
                disp([gap,gapold])
                nonconvind(Qbeta,Qrho,vss,incproc)=mean(nonconv(find(~isnan(nonconv))));
                Coalnonsust(Qbeta,Qrho,vss,incproc)=mean(coalnonsust(find(~isnan(coalnonsust))));
                
                % if (min(min(nonconv))<=0 ||min(min(coalnonsust))<0) & coaldev==1 & vss>1 & itgap>=3
                % if min(min(coalnonsust))<0 & sum(sum(sum(wsbnewsol-wsb)))<0 & coaldev==1 & vss>1 & itgap>=3
                %disp('Coalition not sustainable')
                %itgap
                %sum(sum(sum(wsbnewsol-wsb)))
                % [nonconv,S]
                % pause
                % break
                
                %end
                if  (min(min(nonconv))<=0 ||min(min(coalnonsust))<0) & itgap>=3 & coaldev==1
                    itgap
                    disp('Autarky or coalition not sustainable')
                    [nonconv,S]
                    for j=1:state
                        csol(:,j,1)=caut(j,1);
                        csol(:,j,2)=caut(j,2);
                    end
                    % pause
                    %break
                end
                
                csolold=csol;
                wsb=wsbnewsol;
                ewsb(:,:,1)=(P*wsb(:,:,1)')'; % tessa: great spot!
                ewsb(:,:,2)=(P*wsb(:,:,2)')';
                
                if length(gaphist)>50 & gap>min(gaphist)+0.05
                    gapincrease =1;
                end
            end
            
            CRITGAP(vss)=gaphist(itgap-1);
            
            
            
            if itgap>=maxitgap
                nonconvind(Qbeta,Qrho,vss,incproc)=gap;
            end
            
            if gapincrease==1
                nonconvind(Qbeta,Qrho,vss,incproc)=-gap;
            end
            % =============================================================
            % calculate values of outside option for coalitional deviations
            % =============================================================
            % if insurance group of size vss was unsustainable, use
            % previous value for outside option
            
            if coaldev==1 && (nonconvind(Qbeta,Qrho,vss,incproc)<=0 || Coalnonsust(Qbeta,Qrho,vss,incproc)<=0) && itgap>=3 && vss>1
                % this is redundant, but shows you that if vss is
                % unsustainable, the outside option for vss+1 remains that
                % for vss
                waut_lag(:,1)=waut_lag(:,1);
                waut_lag(:,2)=waut_lag(:,2);
                
                
                
                P_lag=P_lag;
                S_lag=S_lag;
                % if the current insurance group is not sustainable, it
                % cannot be the outside option for the next bigger group,
                % so increase the 'lagind' indicator by one
                lagind=lagind+1;
                gap=-1;
                STABLE(vss)=0;
                % ... or if there is no previous, use just conditional exp of
                % autarky consumption as calculated above
            elseif (nonconvind(Qbeta,Qrho,vss,incproc)<=0 || Coalnonsust(Qbeta,Qrho,vss,incproc)<=0)  && itgap>=3 && (vss==1 || coaldev==0)
                waut_lag(:,1)=eeewaut(:,1);
                waut_lag(:,2)=eeewaut(:,2);
                P_lag=P;
                S_lag=S;
                lagind=1;
                gap=-1;
                
                %%this shouldn't really happen, but for very low discount
                %%factor and rho, it does sometimes not manage it. YOU
                %%NEED TO CHECK THIS AGAIN CAREFULLY!!!!
                S_rov_KEEP(1:state/3,:,vss,Qrho)=S_rov;
                S_KEEP(1:state,:,vss,Qrho)=S;
                PHI_KEEP(:,1:state,:,vss,Qrho)=phisol;
                %EWSB_KEEP(:,1:state,:,vss)=ewsb;
                %EWFIRSTBEST_KEEP(:,1:state,:,vss)=ewfirstbest;
                %WSB_KEEP(:,1:state,:,vss)=wsb;
                STABLE(vss,Qrho)=1;
                IDIO_DIST_ROV_KEEP(1:state/3,1:3,vss,Qrho)=idio_dist;
            else
                lagind=1;
                % when ins group of size vss is sustainable, its value
                % becomes outside option for next bigger group size
                for j=1:state
                    % outside option for next larger insurance group is
                    % value with equal weights vss ...
                    eewaut(j,1)=wsb((n+1)/2,j,1);
                    eewaut(j,2)=wsb((n+1)/2,j,2);
                    
                end
                waut_lag=eewaut;
                S_lag=S;
                P_lag=P;
                
                
                S_rov_KEEP(1:state/3,:,vss,Qrho)=S_rov;
                S_KEEP(1:state,:,vss,Qrho)=S;
                PHI_KEEP(:,1:state,:,vss,Qrho)=phisol;
                %EWSB_KEEP(:,1:state,:,vss)=ewsb;
                %EWFIRSTBEST_KEEP(:,1:state,:,vss)=ewfirstbest;
                %WSB_KEEP(:,1:state,:,vss)=wsb;
                STABLE(vss)=1;
                IDIO_DIST_ROV_KEEP(1:state/3,1:3,vss,Qrho)=idio_dist;
            end
            if itgap>=maxitgap
                nonconvind(Qbeta,Qrho,vss,incproc)=gap;
            end
            % WAUT_KEEP(1:state,1:2,vss)=waut;
            
            %%some diagnostics
            if vss>2
                diagnostics(1,1:2,vsscount,Qrho)=mean(ewfirstbest(451,:,:));
                diagnostics(2,1:2,vsscount,Qrho)=mean(ewsb(451,:,:));
                diagnostics(3,1:2,vsscount,Qrho)=mean(waut);
                rr=find(genuine<0);
                if size(rr)>0
                    diagnostics(4,1,vsscount)=mean(genuine(rr));
                    
                    
                    GENUINE(vsscount,Qbeta,Qrho)=mean(genuine(rr));
                else
                    diagnostics(4,1,vsscount)=0;
                    
                    
                    GENUINE(vsscount,Qbeta,Qrho)=0;
                end
            end
        end
        
    end
    
    Nvill=1;
    
    var_dyagglog=zeros(Nvill*(vss+1),vss_high);var_dyvillagglog=zeros(Nvill*(vss+1),vss_high);vardaggylog=zeros(Nvill*(vss+1),vss_high);betayc=zeros(2,Nvill*(vss+1));vardcshare=zeros(Nvill*(vss+1),vss_high);vardcshare_dylow=zeros(Nvill*(vss+1),vss_high);vardcshare_dyhigh=zeros(Nvill*(vss+1),vss_high);varclog=zeros(Nvill*(vss+1),vss_high);vardclog=zeros(Nvill*(vss+1),vss_high);varylog=zeros(Nvill*(vss+1),vss_high);vardylog=zeros(Nvill*(vss+1),vss_high);varclog_cond=zeros(Nvill*(vss+1),vss_high);varylog_cond=zeros(Nvill*(vss+1),vss_high);
    varc_yhigh_cond=zeros(Nvill*(vss+1),vss_high);varc_ylow_cond=zeros(Nvill*(vss+1),vss_high);vardclog_cond=zeros(Nvill*(vss+1),vss_high);vardylog_cond=zeros(Nvill*(vss+1),vss_high); varc_yhigh=zeros(Nvill*(vss+1),vss_high);
    vardc_dypos_cond=zeros(Nvill*(vss+1),vss_high);
    vardc_dyneg_cond=zeros(Nvill*(vss+1),vss_high); varc_yhigh=zeros(Nvill*(vss+1),vss_high); varc_ylow=zeros(Nvill*(vss+1),vss_high);vardc_dypos=zeros(Nvill*(vss+1),vss_high); vardc_dyneg=zeros(Nvill*(vss+1),vss_high);
    
    betadcdy_agg=zeros(3,Nvill*(vss+1));   betadcdy=zeros(2,vss_high);  betadcdy_high=zeros(1,Nvill*(vss+1));betadcdy_low=zeros(1,Nvill*(vss+1));relcvar=zeros(Nvill*(vss+1),vss_high); reldcvar=zeros(Nvill*(vss+1),vss_high); reldc0var=zeros(Nvill*(vss+1),vss_high); relcvar_cond=zeros(Nvill*(vss+1),vss_high); reldcvar_cond=zeros(Nvill*(vss+1),vss_high);
    beta_cond=zeros(2,Nvill*(vss+1),vss_high); beta_ylow_cond=zeros(4,Nvill*(vss+1),vss_high); beta_yhigh_cond=zeros(4,Nvill*(vss+1),vss_high) ;
    beta_agg=zeros(2,Nvill*(vss+1),vss_high); beta_yhigh_agg=zeros(2,Nvill*(vss+1),vss_high); beta_ylow_agg=zeros(2,Nvill*(vss+1),vss_high);
    IDIO_DIST=[ones(state/3,1)*[1 0 0] + idio_dist; ones(state/3,1)*[0 1 0] + idio_dist; ones(state/3,1)*[0 0 1]+idio_dist];
    
    distribution=size(IDIO_DIST,1);
    
    
    
    
    
    %=============================================
    % Start simulation
    % ============================================
    % number of insurance groups equals village size
    % vss_high+1 divided by size of ins group vss+1
    s = RandStream('mt19937ar','Seed',0)
    RandStream.setGlobalStream(s);
    Npanel=Nvill*(vss+1);
    csim0=nan(Nvill,vss+1,T);
    income0=nan(Nvill,vss+1,T);
    income0ind=nan(Nvill,vss+1,T);
    phisim0=nan(Nvill,vss+1,T);
    Yvill=nan(1,T);
    Yagg=nan(Nvill,T);
    Ssim=nan(Nvill,vss+1,T);
    gammarfun=nan(Nvill,vss+1,T);
    phinew=nan(Nvill,vss+1,T);
    indSjj=nan(T,vss+1);
    
    disp('simulating')
    stream = RandStream.getGlobalStream
    reset(stream);
    shockinittemp=randi(stream,3,Nvill*(vss+1),1);
    shockinit=reshape(shockinittemp,Nvill,vss+1);
    %  shockinit=randi(3,Nvill,vss+1);
    income0ind(:,:,1)=shockinit;
    income0(:,:,1)=incomevals(income0ind(:,:,1));
    csim0(:,:,1)=income0(:,:,1);
    
    %==========================================================
    % Simulate the economy: CD only for the largest group size;
    % ID for all
    % =========================================================
    % This lists all possible combinations of individual income
    % values (IDIO_DIST), and the aggregate states (pair of individual income and ROV-income) associated
    % with every individual income value
    
    stable=(find(STABLE==1));
    
    if coaldev==1                       % for coalitional deviations, all the stable ones
        vss_set_sim=stable;
    else                                % for individual deviations, buildup (but with holes)
        vss_set_sim=vss_set;
        
    end
    
    
    
    for vss=[vss_set_sim]
        
        state=STATE(vss);
        if rov_average==1
            ScaleRov=vss;
        else
            ScaleRov=1;
        end
        
        
        
        
        disp('Reminder: now simulating for village,coaldev,incproc')
        disp([village coaldev incproc])
        disp('now simulating for parameters')
        disp([Qbeta,Qrho])
        disp('vss')
        disp(vss)
        
        
        tic
        
        %initial state
        t=1;
        % initial values of income and consumption (set =
        % income)
        
        
        
        Yagg(:,t)=sum(income0(:,:,t),2);
        
        
        IDIO_DIST_sim=[];
        
        while t<T
            t=t+1;
            %%simulate next period income taking into account
            %% persistence across time
            shock=rand(stream,1,vss+1);
            
            for k=1:vss+1
                
                rr=(find(shock(:,k)<=cumQ(income0ind(:,k,t-1),1)));
                income0ind(rr,k,t)=1;
                rr=(find(cumQ(income0ind(:,k,t-1),1)<shock(:,k) & cumQ(income0ind(:,k,t-1),2)>=shock(:,k)));
                income0ind(rr,k,t)=2;
                rr=(find(cumQ(income0ind(:,k,t-1),2)<shock(:,k) & cumQ(income0ind(:,k,t-1),3)>=shock(:,k)));
                income0ind(rr,k,t)=3;
                income0(:,k,t)=incomevals(income0ind(:,k,t));
            end
            
            Yagg(:,t)=sum(income0(:,:,t),2);
            gammarfun(:,:,t)=(Yagg(:,t-1)*ones(1,vss+1)-csim0(:,:,t-1))./csim0(:,:,t-1)/ScaleRov;
            gammarfun(:,:,t)=min(gammarfun(:,:,t),max(gammagrid(:,4))*ones(Nvill,vss+1));
            
            %Ssim(:,:,t)=indS(IDIO_DIST_sim(:,1),:);
            
            
            IDIO_DIST_sim=find(((IDIO_DIST(:,1,1)==(sum(income0ind(1,:,t)==1))).*(IDIO_DIST(:,2,1)==(sum(income0ind(1,:,t)==2))).*(IDIO_DIST(:,3,1)==(sum(income0ind(1,:,t)==3))))==1)';
            for k=1:vss+1
                % choose for individual k the aggregate state
                % that applies
                %indSjj(t,k)=IDIO_DIST_sim(S(IDIO_DIST_sim,1)<=income0(ivill,k,t)+eps & S(IDIO_DIST_sim,1)>=income0(ivill,k,t)-eps);
                indSjj(t,k)=IDIO_DIST_sim(S(IDIO_DIST_sim,1)==income0(1,k,t));
                % check that aggregate and individual income0s are
                % consistent with the states caculated above
                if abs(Yagg(1,t)-(S(indSjj(t,k),1)+ScaleRov*S(indSjj(t,k),2)))>0.000001
                    stop
                end
                %%and find the gammarfun2
                
            end
            Ssim(1,:,t)=indSjj(t,:);
            
            countrho=0;
            for Qrho=1:size(rhovec,2)
                rho=rhovec(Qrho);
                
                S=S_KEEP(1:state,:,vss,Qrho);
                S_rov=S_rov_KEEP(1:state/3,:,vss,Qrho);
                idio_dist=IDIO_DIST_ROV_KEEP(1:state/3,1:3,vss,Qrho);
                phisol=PHI_KEEP(:,1:state,:,vss,Qrho);
                for k=1:(vss+1)/length(rhovec)
                    countrho=countrho+1;
                    BLA=interp1(gammagrid(:,4),phisol(:,Ssim(:,countrho,t),1),gammarfun(:,countrho,t));
                    phinew(:,countrho,t)=diag(BLA);
                    rho_longvec(countrho)=rhovec(Qrho);
                    %PHI_temp(countrho)=(1+phinew(:,1,t)*ones(1,vss+1))).^(1/rhovec(Qrho));
                    margut_lag(:,countrho)=csim0(:,countrho,t-1)^(-rhovec(Qrho));
                end
            end
            x0=csim0(:,:,t-1);
            [csim0(:,:,t),output]=conssolve_het(x0,phinew(:,:,t),Yagg(:,t),margut_lag,rho_longvec);
            % group=[1:vss+1];
            %for ivill=1:Nvill
            %   phinewother(ivill,1:vss+1,t)=0;
            %for k=1:vss+1
            %     gammarfun2(ivill,:,t)=csim0(ivill,:,t-1)./csim0(ivill,k,t-1);
            
            %   numerator=0;
            %  denominator=0;
            % rr=(group~=k);
            %crov=0;
            %for gg=group(rr)
            %  numerator=numerator + (1+phinew(ivill,gg,t))^(1/rho)*gammarfun2(ivill,gg,t);
            % denominator=denominator + gammarfun2(ivill,gg,t);
            %crov=crov+csim0(ivill,gg,t-1);
            %end
            %phinewother(ivill,k,t)= (numerator/denominator)^(rho) -1;
            %Crov(ivill,k,t-1)=crov;
            
            %end
            %end
            
            %if (nonconvind(Qbeta,Qrho,vss,incproc)<=0 || Coalnonsust(Qbeta,Qrho,vss,incproc)<=0)
            %%now calculate the gammaraux for n individuals,
            %%with 1 as the reference individual in order to
            %%solve for consumption directly
            %csim0(:,:,t)=income0(:,:,t);
            %else
            %(margut_lag(:,1)./margut_lag(:,2:end))=(1+phinew(:,2:end,t)./(1+phinew(:,1,t).*(margut_lag(:,1)./margut_lag(:,2:end));
            %gammaraux=((1+phinew(:,:,t))./(1+phinew(:,1,t)*ones(1,vss+1))).^(1/rho).*csim0(:,:,t-1)./(csim0(:,1,t-1)*ones(1,vss+1));
            
            %gammaraux=((1+phinew(:,:,t))./(1+phinew(:,1,t)*ones(1,vss+1))).^(1/rho).*csim0(:,:,t-1)./(csim0(:,1,t-1)*ones(1,vss+1));
            %%and new consumption
            %for ivill=1:Nvill
            %for k=1:vss+1
            %csim0_2(ivill,k,t)=Yagg(ivill,t)/(1+((1+phinewother(ivill,k,t))/(1+phinew(ivill,k,t))).^(1/rho)*Crov(ivill,k,t-1)/csim0_2(ivill,k,t-1));
            %end
            %end
            
            % csim0(:,:,t)=(gammaraux./(sum(gammaraux,2)*ones(1,vss+1))).*(Yagg(:,t)*ones(1,vss+1));
            %end
        end
    end
    disp('simulating takes')
    tsim(Qbeta,Qrho)=toc
    
    tic
    %reshape
    
    AA=[];
    BB=[];
    for t=1:T
        AA=[AA reshape(csim0(:,:,t)',Nvill*(vss+1),1)];
        BB=[BB reshape(income0(:,:,t)',Nvill*(vss+1),1)];
    end
    csim0=AA;
    
    income0=BB;
    
    disp('reshaping takes')
    toc
    
    
    
    % =========================================================
    % calculate moments of consumption and income growth
    % =========================================================
    csim =[]; cshare_raw =[]; dcshare_raw =[]; cc =[]; ccc =[]; yy =[]; aggyy =[]; aggvillincome=[]; yyy_cond  =[]; dc_log_cond =[]; dy_log_cond =[]; ccc_cond=[];
    cerr=randn(size(csim0)); logcsim0_groupmeandev=[]; logcsim_groupmeandev=[]; logincome_groupmeandev=[];
    incerr=randn(size(csim0));   Momentsest1=nan(27,11); ind=[];
    % loop over size of cons measurement error
    for kkk=1:size(csim0,1)
        tempvar(kkk)=var(log(csim0(kkk,:)));
    end
    Tempvar=mean(tempvar);
    for ic=1:N_cerr+1
        csim=exp(log(csim0)+((ic-1)/N_cerr*maxcerr/(1-(ic-1)/N_cerr*maxcerr)*Tempvar)^0.5*cerr);
        logcsim=log(csim);
        logcsim0=log(csim0);
        logcsim0_meandev=logcsim0-ones(size(csim0,1),1)*mean(logcsim0,1);
        logcsim_meandev=logcsim-ones(size(csim0,1),1)*mean(logcsim,1);
        % add inc meas error back in
        income=exp(log(income0)+varnu(incproc)^0.5*incerr);
        logincome=log(income);
        logincome0=log(income0);
        logincome_meandev=logincome-ones(size(csim0,1),1)*mean(logincome,1);
        income0_meandev=income0-ones(size(csim0,1),1)*mean(income0,1);
        
        
        for Qrho=1:length(rhovec)
            ngroup=(vss+1)/length(rhovec);
            aggvillincome((Qrho-1)*(ngroup)+1:Qrho*(ngroup),1:T)=ones(ngroup,1)*sum(income((Qrho-1)*(ngroup)+1:Qrho*(ngroup),:),1);
            logaggvillincome=log(aggvillincome);
            logcsim0_groupmeandev((Qrho-1)*(ngroup)+1:Qrho*(ngroup),1:T)=logcsim0((Qrho-1)*(ngroup)+1:Qrho*(ngroup),:)-ones(ngroup,1)*mean(logcsim0((Qrho-1)*(ngroup)+1:Qrho*(ngroup),:),1);
            logcsim_groupmeandev((Qrho-1)*(ngroup)+1:Qrho*(ngroup),1:T)=logcsim((Qrho-1)*(ngroup)+1:Qrho*(ngroup),:)-ones(ngroup,1)*mean(logcsim((Qrho-1)*(ngroup)+1:Qrho*(ngroup),:),1);
            logincome_groupmeandev((Qrho-1)*(ngroup)+1:Qrho*(ngroup),1:T)=logincome((Qrho-1)*(ngroup)+1:Qrho*(ngroup),1:T)-ones(ngroup,1)*mean(logincome((Qrho-1)*(ngroup)+1:Qrho*(ngroup),:),1);
        end
        aggincome=sum(income,1);
        
        
        
        % === discard first Tinit-1 periods
        % this chooses log or level consumption for
        % moments
        if momentclevinlog==1
            ccc=logcsim(:,Tinit:T);
            ccc_cond=logcsim_meandev(:,Tinit:T);
            ccc_groupcond=logcsim_groupmeandev(:,Tinit:T);
        else
            ccc=(csim(:,Tinit:T));
            ccc_cond=(csim_meandev(:,Tinit:T));
        end
        yyy=logincome(:,Tinit:T);
        yyy_cond=logincome_meandev(:,Tinit:T);
        yyy_groupcond=logincome_groupmeandev(:,Tinit:T);
        
        cc=logcsim(:,Tinit:T)-logcsim(:,Tinit-1:T-1);
        cc_cond=logcsim_meandev(:,Tinit:T)-logcsim_meandev(:,Tinit-1:T-1);
        cc_groupcond=logcsim_groupmeandev(:,Tinit:T)-logcsim_groupmeandev(:,Tinit-1:T-1);
        
        yy=logincome(:,Tinit:T)-logincome(:,Tinit-1:T-1);
        yy_cond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
        yy_groupcond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
        aggyy=log(sum(income(:,Tinit:T),1))-log(sum(income(:,Tinit-1:T-1),1));
        aggvillyy=log(aggvillincome(:,Tinit:T))-log(aggvillincome(:,Tinit-1:T-1));
        
        %for ivill=1:Nvill
        %aggvillyy=log(aggvillincome(ivill,Tinit:T))-log(aggvillincome(ivill,Tinit-1:T-1));
        %var_dyvillagglog(ivill,vss)=var(aggvillyy(ivill,:));
        %end
        vardaggylog(1,vss)=var(log(aggincome(1,Tinit:T))-log(aggincome(1,Tinit-1:T-1)));
        
        % calculate moments for each individiual
        for ind=1:size(csim,1)
            
            dc_log_cond(ind,:)=(ccc_cond(ind,2:end))-(ccc_cond(ind,1:end-1));
            dy_log_cond(ind,:)=(yyy_cond(ind,2:end))-(yyy_cond(ind,1:end-1));
            dc_log_groupcond(ind,:)=(ccc_groupcond(ind,2:end))-(ccc_groupcond(ind,1:end-1));
            dy_log_groupcond(ind,:)=(yyy_groupcond(ind,2:end))-(yyy_groupcond(ind,1:end-1));
            
            % conditioning on income increases, and income
            % decreases
            % NB: we allocate no change to decrease!!
            rrhc =[]; rrlc =[]; rrl =[]; rrh=[];
            rrhc=find(yy(ind,:)>0);
            rrlc=find(yy(ind,:)<=0);
            rrhc_cond=find(yy_cond(ind,:)>0);
            rrlc_cond=find(yy_cond(ind,:)<=0);
            rrhc_groupcond=find(yy_groupcond(ind,:)>0);
            rrlc_groupcond=find(yy_groupcond(ind,:)<=0);
            
            
            %%THE VARIANCES
            
            vardclog(ind,vss)=var(cc(ind,:));
            vardylog(ind,vss)=var(yy(ind,:));
            vardc_dypos(ind,vss)=var(cc(ind,rrhc));
            vardc_dyneg(ind,vss)=var(cc(ind,rrlc));
            
            %%conditional on village income!
            vardclog_cond(ind,vss)=var(dc_log_cond(ind,:));
            vardylog_cond(ind,vss)=var(dy_log_cond(ind,:));
            vardc_dypos_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)>0));
            vardc_dyneg_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)<=0));
            
            %%conditional on group income!
            %Step 1: create your residuals
            vardclog_groupcond(ind,vss)=var(dc_log_groupcond(ind,:));
            vardylog_groupcond(ind,vss)=var(dy_log_groupcond(ind,:));
            vardc_dypos_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)>0));
            vardc_dyneg_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)<=0));
            vardc_dypos(ind,vss)=var(cc(ind,rrhc));
            vardc_dyneg(ind,vss)=var(cc(ind,rrlc));
            
            %%and the relative variances
            
            reldcvar(ind,vss)=var(cc(ind,rrhc))-var(cc(ind,rrlc));
            
            reldcvar_cond(ind,vss)=vardc_dypos_cond(ind,vss)-vardc_dyneg_cond(ind,vss);
            reldcvar_groupcond(ind,vss)=vardc_dypos_groupcond(ind,vss)-vardc_dyneg_groupcond(ind,vss);
            
            
            %%THE BETA COEFFICIENTS
            
            %1) Unconditional
            %risk-sharing whole sample
            covtemp=cov(cc(ind,:),aggyy(1,:));
            betadcdy_agg(ind,vss)=covtemp(2,1)/var(aggyy(1,:));
            covtemp=cov(cc(ind,:),yy(ind,:));
            betadcdy(ind,vss)=covtemp(2,1)/(var(yy(ind,:)));
            %risk-sharing for those with income gains versus
            %those with income losses
            YY=cc(ind,:)';
            yy_hc(ind,:)=zeros(1,size(cc,2));
            yy_hc(ind,rrhc)=yy(ind,rrhc);
            oneone=ones(1,size(cc,2));
            ones_hc(ind,:)=zeros(1,size(cc,2));
            ones_hc(ind,rrhc)=oneone(rrhc);
            XX=[ones(1,size(cc,2))',ones_hc(ind,:)',yy(ind,:)',yy_hc(ind,:)'];
            beta_yhigh(:,ind,vss)=regress(YY,XX);
            betadcdy_low(ind,vss)=beta_yhigh(3,ind,vss);
            betadcdy_high(ind,vss)=beta_yhigh(4,ind,vss)+beta_yhigh(3,ind,vss);
            
            %2) conditional on villag income
            covtemp=[];
            X=[ones(size(yy_cond(ind,:),2),1),yy_cond(ind,:)'];%
            YY=cc_cond(ind,:)';
            beta_cond(:,ind,vss)=inv(X'*X)*(X'*YY);
            
            X=[ones(size(aggyy(1,:),2),1),aggyy'];%
            YY=cc(ind,:)';
            beta_agg(:,ind,vss)=inv(X'*X)*(X'*YY);
            
            %risk-sharing for those with income gains versus
            %those with income losses
            YY_cond=cc_cond(ind,:)';
            yy_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
            yy_hc_cond(ind,rrhc_cond)=yy_cond(ind,rrhc_cond);
            oneone=ones(1,size(cc_cond,2));
            ones_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
            ones_hc_cond(ind,rrhc_cond)=oneone(rrhc_cond);
            XX_cond=[ones(1,size(cc_cond,2))',ones_hc_cond(ind,:)',yy_cond(ind,:)',yy_hc_cond(ind,:)'];
            beta_yhigh_cond(:,ind,vss)=regress(YY_cond,XX_cond);
            
            %
            
            %2) conditional on group income
            covtemp=[];
            X=[ones(size(yy_groupcond(ind,:),2),1),yy_groupcond(ind,:)'];%
            YY=cc_groupcond(ind,:)';
            beta_groupcond(:,ind,vss)=inv(X'*X)*(X'*YY);
            
            X=[ones(size(aggvillyy(ind,:),2),1),aggvillyy(ind,:)'];%
            YY=cc(ind,:)';
            beta_groupagg(:,ind,vss)=inv(X'*X)*(X'*YY);
            
            %risk-sharing for those with income gains versus
            %those with income losses
            YY_groupcond=cc_groupcond(ind,:)';
            yy_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
            yy_hc_groupcond(ind,rrhc_groupcond)=yy_groupcond(ind,rrhc_groupcond);
            oneone=ones(1,size(cc_groupcond,2));
            ones_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
            ones_hc_groupcond(ind,rrhc_groupcond)=oneone(rrhc_groupcond);
            XX_groupcond=[ones(1,size(cc_groupcond,2))',ones_hc_groupcond(ind,:)',yy_groupcond(ind,:)',yy_hc_groupcond(ind,:)'];
            beta_yhigh_groupcond(:,ind,vss)=regress(YY_groupcond,XX_groupcond);
            
            
            
            
            
            
            
            
            
        end
        var_dyagglog(1,vss)=var(aggyy);
        % vector of moments
        % col 1-5 are parameters
        Momentsest1(1:5,ic)=[vss,beta,rho,rho1vec(incproc),varnu(incproc)]';
        % 6-11 relative variance dc/dy, dc(dy>0)/dc(dy<0);  unconditional (6:7), conditional on village income
        % (8:9), conditonal on group income (10:11)
        Momentsest1(6:11,ic)=[mean(vardclog(1:size(csim,1),vss))/mean(vardylog(1:size(csim,1),vss)), mean(reldcvar(1:size(csim,1),vss))...
            mean(vardclog_cond(1:size(csim,1),vss))/mean(vardylog_cond(1:size(csim,1),vss)), mean(reldcvar_cond(1:size(csim,1),vss))...
            mean(vardclog_groupcond(1:size(csim,1),vss))/mean(vardylog_groupcond(1:size(csim,1),vss)), mean(reldcvar_groupcond(1:size(csim,1),vss))]';
        % 12-14: beta, beta high, beta low unconditional
        Momentsest1(12:14,ic)=[mean(betadcdy(1:size(csim,1),vss)),mean(betadcdy_high(1:size(csim,1),vss)),mean(betadcdy_low(1:size(csim,1),vss))]';
        % 13_16: beta, beta high, beta low conditional on
        % village income
        Momentsest1(15:17,ic)=[mean(beta_cond(2,1:size(csim,1),vss),2),mean(beta_yhigh_cond(3,1:size(csim,1),vss),2)+mean(beta_yhigh_cond(4,1:size(csim,1),vss),2),mean(beta_yhigh_cond(3,1:size(csim,1),vss),2)]';
        % 18:20: beta, beta high, beta low conditional on
        % group income
        Momentsest1(18:20,ic)=[mean(beta_groupcond(2,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(3,1:size(csim,1),vss),2)+mean(beta_yhigh_groupcond(4,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(3,1:size(csim,1),vss),2)]';
        %21:22 are variances of dclog, dylog
        
        Momentsest1(21:22,ic)=[mean(vardclog(1:size(csim,1),vss)),mean(vardylog(1:size(csim,1),vss))];
        %23:28, relative variance dc/dy (dy>0) and var dc/dy(dy<0);
        %unconditinoal 23:24, conditional on village income 25:26,
        %conditional on group income 27:28
        Momentsest1(23:24,ic)=[mean(vardc_dypos(1:size(csim,1),vss));mean( vardc_dyneg(1:size(csim,1),vss))];
        Momentsest1(25:26,ic)=[mean(vardc_dypos_cond(1:size(csim,1),vss));mean(  vardc_dyneg_cond(1:size(csim,1),vss))];
        Momentsest1(27:28,ic)=[mean(vardc_dypos_groupcond(1:size(csim,1),vss));mean( vardc_dyneg_groupcond(1:size(csim,1),vss))];
        
        %%and the critical gap
        
        Momentsest1(29,1)=CRITGAP(vss);
        Momentsest1(30,ic)=[mean( vardylog_cond(1:size(csim,1),vss))];
        Momentsest1(31,ic)=[rho1vec(incproc)];
        countrho=0;
        for Qrho=1:length(rhovec)
            countrho=countrho+1;
            ngroup=(vss+1)/length(rhovec);
            Momentsest1(31+(countrho-1)*5+1:31+(countrho)*5,ic)=[mean(vardclog_cond((Qrho-1)*(ngroup)+1:Qrho*(ngroup),vss))/mean(vardylog_cond((Qrho-1)*(ngroup)+1:Qrho*(ngroup),vss))
                mean(beta_cond(2,(Qrho-1)*(ngroup)+1:Qrho*(ngroup),vss))
                mean(reldcvar_cond((Qrho-1)*(ngroup)+1:Qrho*(ngroup),vss))
                mean(beta_yhigh_cond(3,(Qrho-1)*(ngroup)+1:Qrho*(ngroup),vss),2)+mean(beta_yhigh_cond(4,(Qrho-1)*(ngroup)+1:Qrho*(ngroup),vss),2)
                mean(beta_yhigh_cond(3,(Qrho-1)*(ngroup)+1:Qrho*(ngroup),vss),2)]
            
        end
    end
    
    
    
    
    toc
    
    
    if coaldev==0
        Momentsest(:,:,Qbeta,1,vss)=[Momentsest1];
        %                Income0(:,:,Qbeta,Qrho,vss)=income0(1:vss_high_inddev,:);
        %  Csim0(:,:,Qbeta,Qrho,vss)=csim0(1:vss_high_inddev,:);
    elseif coaldev==1
        Momentsest(:,:,Qbeta,1,vss)=[Momentsest1];
        %  Income0(:,:,Qbeta,Qrho,vss)=income0(1:vss_high_inddev,:);
        %Csim0(:,:,Qbeta,Qrho,vss)=csim0(1:vss_high_inddev,:);
    end
    % pause
    if buildup==0
        
        cd(outputpath)
        save(filetosaveest,'Momentsest');
        cd(programpath)
    elseif buildup==1
        save('Momentsestbuildup','Momentsest');
    end
end


cd(outputpath)
eval(['save Het', filetosave, int2str(village),int2str(coaldev), int2str(incproc),int2str(unconditional) '.mat'])
cd(programpath)

%forensics


